// "inflate.cpp" - Implements the inflate 2d stroke procedure given by
// T. Igarashi et al. in the 1999 Siggraph paper:
// "Teddy: A Sketching Interface for 3D Freeform Design"
//
// May be used by itself with only triangle.c/triangle.h, as it has
// no other dependences other than STL.
//
// The code in its current state is a little messy, but works.
// Will clean it up provided someone has a need for it.
/* OpenTeddy 3d sketch modeler based on the paper by T. Igarashi et al.
*      "Teddy: A Sketching Interface for 3D Freeform Design"
*
* Copyright (c) 2004 Alex Ni
*
*
* Permission is hereby granted, free of charge, to any person obtaining
* a copy of this software and associated documentation files (the
* "Software"), to deal in the Software without restriction, including
* without limitation the rights to use, copy, modify, merge, publish,
* distribute, sublicense, and/or sell copies of the Software, and to
* permit persons to whom the Software is furnished to do so, subject to
* the following conditions:
* 
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
* 
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
* MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
* IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
* CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
* TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
* SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/


// C includes
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
extern "C"
{
#include "triangle.h"
}

// C++ includes
#include <list>
#include <vector>
using std::list;
using std::vector;
using std::less;
using std::unary_function;
using std::pair;

#include "otconfig.h"
#include "inflate.h"

void ot_inflate::do_triangulate(double* i_points, size_t i_ptct)
{
	OT_ASSERT(!triangulated);
	OT_ASSERT(i_ptct > 0);
	points = new double[2 * (i_ptct + 1)];  // must assure an extra one in case we have no tri w/ 3 neighbors
	for (int i = 0;i < (2*(i_ptct + 1));i++)
	{
		points[i] = (i < (2 * i_ptct)) ? i_points[i] : 0.0;
	}
	ptct = i_ptct;
	
	struct triangulateio in;
	memset((void*)&in, 0, sizeof(struct triangulateio));
	if (ptct < 3) return ;
	in.numberofpoints = ptct;
	in.numberofpointattributes = 0;
	in.pointlist = (REAL *) points;
	in.pointmarkerlist = NULL;
	
	in.numberofsegments = ptct;
	in.segmentlist = new int[2 * (in.numberofsegments)];
	for (int i = 0;i < in.numberofsegments;i++)
	{
		in.segmentlist[i*2] = i;
		in.segmentlist[(i*2) + 1] = (i + 1) % (in.numberofsegments);
	}
	/* Needed only if segments are output (-p or -c) and -P and -B not used: */
	
	
	in.numberofholes = 0;
	in.numberofregions = 0;
	in.regionlist = NULL; //(REAL *) malloc(in.numberofregions * 4 * sizeof(REAL));
	//  in.regionlist[0] = 0.5;
	//  in.regionlist[1] = 5.0;
	//  in.regionlist[2] = 7.0;            /* Regional attribute (for whole mesh). */
	//  in.regionlist[3] = 0.1;          /* Area constraint that will not be used. */
	// i have no idea what these mean, might want to use them later - alexni
	
	memset((void*)&out, 0, sizeof(struct triangulateio));
	
	out.pointlist = (REAL *) NULL;            /* Not needed if -N switch used. */
	/* Not needed if -N switch used or number of point attributes is zero: */
	out.pointattributelist = (REAL *) NULL;
	out.pointmarkerlist = (int *) NULL; /* Not needed if -N or -B switch used. */
	out.trianglelist = (int *) NULL;          /* Not needed if -E switch used. */
	/* Not needed if -E switch used or number of triangle attributes is zero: */
	out.triangleattributelist = (REAL *) NULL;
	out.neighborlist = (int *) NULL;         /* Needed only if -n switch used. */
	/* Needed only if segments are output (-p or -c) and -P not used: */
	out.segmentlist = (int *) NULL;
	/* Needed only if segments are output (-p or -c) and -P and -B not used: */
	out.segmentmarkerlist = (int *) NULL;
	out.edgelist = (int *) NULL;             /* Needed only if -e switch used. */
	out.edgemarkerlist = (int *) NULL;   /* Needed if -e used and -B not used. */
	triangulated = true;
	triangulate("pzeNn", &in, &out, (struct triangulateio*) NULL);
	// p = input is segments, not just points (so we can have a concave polygon)
	// z = start index at 0 (??? what is the alternative?)
	// e = output edges.
	// N = don't copy in.pointlist to out.pointlist (since it is this same list w/ the other flags)
	// n = neighbors are listed (saves us calculating some connectivity info)
	int j = 1;
	
	delete[] in.segmentlist;
	
	//fields in out we eventually need to get rid of:
	/*
	free(out.trianglelist);
	free(out.neighborlist);
	free(out.segmentlist);
	free(out.segmentmarkerlist);
	free(out.edgelist);
	free(out.edgemarkerlist);
	*/
}


int ot_inflate::add_npoint(double p1, double p2, bool external)
{
	if (external)
	{
		npoints[ptct*2] = p1;
		npoints[ptct*2 + 1] = p2;
		return ptct++;
	}
	
	npoints[nptct*2] = p1;
	npoints[nptct*2 + 1] = p2;
	return nptct++;
}

#define NEIGH(i,j)   out.neighborlist[ ((i)*3) + (j) ]
#define PT(i,j)      points[ (2*(i)) + (j) ]
#define NPT(i,j)     npoints[ (2*(i)) + (j) ]
#define N3PT(i,j)	 n3pts[ (3*(i)) + (j) ]
#define EDGE(i,j)    out.edgelist[ ((i)*2) + (j) ]
#define TRI(i,j)     out.trianglelist[ (3*(i)) + (j) ]
#define NEIGHCNT(a)  (  ( (NEIGH((a),0) == -1) ? 0 : 1 )+\
	( (NEIGH((a),1) == -1) ? 0 : 1 )+\
( (NEIGH((a),2) == -1) ? 0 : 1 )  )



// shared edge of two triangles
bool ot_inflate::shedge(int t1, int t2, exp_edge* ee )
{
	// if a triangle is -1, they definitely don't share an edge
	if ( (t1 < 0) || (t2 < 0) ) return false;
	
	// for each face of each, for each face of each, try to match up
	for (int i = 0;i < 3;i++)
	{
		for (int j = 0;j < 3;j++)
		{
			if ( ( ( TRI(t1, i) == TRI(t2, j) ) &&
				( TRI(t1, (i + 1) % 3) == TRI(t2, (j + 1) % 3) ) ) ||
				( ( TRI(t1, (i + 1) % 3) == TRI(t2, j) ) &&
				( TRI(t1, i) == TRI(t2, (j + 1) % 3) ) ) )
			{
				ee->p1 = TRI(t1, i);
				ee->p2 = TRI(t1, (i + 1) % 3);
				return true;
			}
		}
	}
	return false;
}

#define SWAP(x,y,t) {t = x; x = y; y = t;}

// return edge and adj triangle of a sleeve triangle with a neighbor (other than ee)
bool ot_inflate::adj_sleeve_edge_tri(int tri, exp_edge* ee, exp_edge* adjee, int* adjt)
{
	for (int i = 0;i < 3;i++)
	{
		if (shedge(tri, *adjt = NEIGH(tri, i), adjee ) && !(*adjee == *ee))
		{
			// make it so ee[0] == *adjee[0]
			if ( (*ee)[0] != (*adjee)[0])
			{
				int temp;
				bool s1, s2;
				s1 = s2 = false;
				
				if ((*ee)[0] == (*adjee)[1]) s1 = true;
				else if ((*ee)[1] == (*adjee)[0]) s2 = true;
				else s1 = s2 = true;
				
				if (s1) SWAP((*adjee)[0], (*adjee)[1], temp)
					if (s2) SWAP((*ee)[0], (*ee)[1], temp)
			}
			return true;
		}
	}
	return false;
}


int ot_inflate::find_root_joint(int tri)
{
	// sanity check
	int k = NEIGHCNT(tri);
	if (tri >= out.numberoftriangles) return -1;
	
	// if this is a joint, we're done
	printf("Tri %d has neih %d\n", tri, NEIGHCNT(tri));
	
	if (NEIGHCNT(tri) == 3)
		return tri;
	// otherwise mark it, try a neighbor
	fjmarked[tri] = 1;
	int foundtri = -1;
	for (int i = 0;i < 3;i++)
	{
		int neigh = NEIGH(tri, i);
		if ( (foundtri == -1) && (neigh != -1) && !fjmarked[neigh] )
		{
			foundtri = find_root_joint(neigh);
		}
	}
	// unmark
	fjmarked[tri] = 0;
	
	return foundtri;
}

static double dist(double x1, double y1, double x2, double y2)
{
	double dx = (x1 - x2);
	double dy = (y1 - y2);
	return sqrt(dx*dx + dy*dy);
}

void ot_inflate::add_edge_to_fanpts(list<int>* fanpts, exp_edge ee)
{
	std::list<int>::iterator begi = fanpts->begin();
	std::list<int>::iterator endi = fanpts->end();
	endi--;
	
	// add to the opposite we find it on
	if (ee[0] == (*begi))
	{
		fanpts->push_back(ee[1]);
	}
	else if (ee[1] == (*begi))
	{
		fanpts->push_back(ee[0]);
	}
	else if (ee[0] == (*endi))
	{
		fanpts->insert(begi, ee[1]);
	}
	else if (ee[1] == (*endi))
	{
		fanpts->insert(begi, ee[0]);
	}
	else
	{
		OT_ASSERT(0);
	}
}

// external points are the ones that are input (the edge of the polygon)
// internal points are on the spine
bool ot_inflate::is_external(int pt)
{
	if (pt < ptct) return true;
	return false;
}


// add edge going from external point to spine
void ot_inflate::add_new_intext_edge(exp_edge ee)
{
	// check if not already added
	if (intext_h.find(ee) != intext_h.end()) return ;
	
	intext_v.push_back(ee);
	intext_h[ee] = intext_v.size() - 1;
	spine_sum_dist[ee[1]] += dist( NPT(ee[0], 0), NPT(ee[0], 1), NPT(ee[1], 0), NPT(ee[1], 1) );
	spine_sum_total[ee[1]] += 1;
}

double ot_inflate::newtri_area(exp_tri etri)
{
	double x0, y0, x1, y1;
	x0 = NPT(etri[1], 0) - NPT(etri[0], 0);
	y0 = NPT(etri[1], 1) - NPT(etri[0], 1);
	x1 = NPT(etri[2], 0) - NPT(etri[0], 0);
	y1 = NPT(etri[2], 1) - NPT(etri[0], 1);
	return (x0*y1) - (y0*x1);
}

void ot_inflate::add_newtri(exp_tri etri)
{
	for (int i = 0;i < 3;i++)
	{
		if (is_external(etri[i]) && !is_external(etri[(i + 1) % 3]))
		{
			add_new_intext_edge( exp_edge( etri[i], etri[(i + 1) % 3] ) );
		}
		if (is_external(etri[(i + 1) % 3]) && !is_external(etri[i]))
		{
			add_new_intext_edge( exp_edge( etri[(i + 1) % 3], etri[i] ) );
		}
	}
	
	if (newtri_area(etri) > 0.0)
		newtri.push_back(etri);
	else
		newtri.push_back( exp_tri( etri[2], etri[1], etri[0] ) );
}

void ot_inflate::build_tri_fin(list<int>* fanpts, int cenpt)
{
	for (std::list<int>::iterator i = fanpts->begin();
	i != fanpts->end(); ++i)
	{
		
		std::list<int>::iterator nextit = i;
		nextit++;
		// skip last (since it's already used in a tri)
		if (nextit == fanpts->end()) continue;
		add_newtri(exp_tri(cenpt, *i, *nextit));
	}
}

// squared weight- tends to fix the problem where one side is longer than the other,
// making the skeleton take a sharp jump in some cases

static void weighted_center(double x1, double y1,
							double x2, double y2,
							double x3, double y3,
							double &cx, double &cy)
{
	double side1, side2, side3;
	double s1x, s1y, s2x, s2y, s3x, s3y;
	side1 = dist(x1, y1, x2, y2);
	side1 = side1 * side1;
	s1x = (x1 + x2) / 2; s1y = (y1 + y2) / 2;
	side2 = dist(x2, y2, x3, y3);
	side2 = side2 * side2;
	s2x = (x2 + x3) / 2; s2y = (y2 + y3) / 2;
	side3 = dist(x3, y3, x1, y1);
	side3 = side3 * side3;
	s3x = (x3 + x1) / 2; s3y = (y3 + y1) / 2;
	
	cx = (s1x * side1 + s2x * side2 + s3x * side3) / (side1 + side2 + side3);
	cy = (s1y * side1 + s2y * side2 + s3y * side3) / (side1 + side2 + side3);
	/*	cx =  (  x1*( abs(x1-x2) + abs(x1-x3) ) + x2*( abs(x2-x1) + abs(x2-x3) )
				+x3*( abs(x3-x1) + abs(x3-x2) ) ) /
				(2*abs(x1-x2) + 2*abs(x1-x3) + 2*abs(x2-x3));
				cy =  (  y1*( abs(y1-y2) + abs(y1-y3) ) + y2*( abs(y2-y1) + abs(y2-y3) )
				+y3*( abs(y3-y1) + abs(y3-y2) ) ) /
				(2*abs(y1-y2) + 2*abs(y1-y3) + 2*abs(y2-y3));
	*/
}


// gen spine triangles
// 'tri' is the tri to generate spine triangles at (old mesh)
// 'ee' is the shared edge between the calling triangle and the cur one (old mesh)

// return points that should be part of the current triangle fan tip (if any)
// or return new (new mesh) edges generated that are replacing 'ee'
//   if sleeve between two joints, that should be null, and
//   those sleeve triangles should just be tesselated

// 'ee' is a struct... we're doing that instead of an edge index since
// otherwise we'd have to keep looking up edge indice from the neighbor list.
gst_return ot_inflate::gen_spine_tris(int tri, exp_edge ee, int force_neigh)
{
	gst_return returned, toret;
	switch (force_neigh ? force_neigh : NEIGHCNT(tri))
	{
	case 1 :    // tip
		{
			// 1. no neighbors to call...
			
			// 2. add tip points to triangle fan list
			// return (0, triangle-fan-list)
			toret.fanpts = new list<int>;
			// make sure to put them in order, center goes in the middle
			int center;
			for (int i = 0;i < 3;i++)
			{
				if ( (TRI(tri, i) != ee[0]) && (TRI(tri, i) != ee[1]) )
					center = TRI(tri, i);
			}
			toret.fanpts->push_back(ee[0]);
			toret.fanpts->push_back(center);
			toret.fanpts->push_back(ee[1]);
			
			toret.newp = -1;
		}
		break;
	case 2:   // sleeve
		{
			
			// 1. call tri on the other side of the sleeve first
			int adjt;
			exp_edge adjee;
			if (!adj_sleeve_edge_tri(tri, &ee, &adjee, &adjt))
			{
				// Uh oh...
				OT_ASSERT(0);
			}
			
			// we're going to need to tesselate this no matter what
			double nx, ny;
			nx = (NPT(adjee[0], 0) + NPT(adjee[1], 0)) / 2;
			ny = (NPT(adjee[0], 1) + NPT(adjee[1], 1)) / 2;
			int tessptadj, tesspt;
			// tessptadj = (nx,ny),
			// allocated as a new point only if we make a triangle fan
			// otherwise it is given to use as a result
			
			tessptadj = tesspt = -1;
			
			returned = gen_spine_tris(adjt, adjee);
			// 2. if triangle fan list is nonzero
			if (returned.fanpts)
			{
				// try circle criterion (point_outside) on called edge 'e',
				// (true if any points lie outside the radius)
				
				// find new point at center of edge (avg of two points)
				double dr;
				dr = dist(nx, ny, NPT(adjee[0], 0), NPT(adjee[0], 1));
				
				bool point_outside = false;
				for (std::list<int>::iterator i = returned.fanpts->begin();
				i != returned.fanpts->end(); ++i)
				{
					if ( dist(nx, ny, NPT(*i, 0), NPT(*i, 1)) > dr )
					{
						point_outside = true;
					}
				}
				
				// if no points outside, we can keep adding to the list.
				if (!point_outside)
				{
					// add new point from the sleeve to the triangle fan list
					toret.fanpts = returned.fanpts;
					// find out which point is new
					add_edge_to_fanpts(toret.fanpts, ee);
					
					// return (0, triangle_Fan_list)
					toret.newp = -1;
					return toret;
					// else
				}
				else
				{
					// create triangle fan at edge
					// add new point at center of edge (avg of two points)
					tessptadj =
						add_npoint( nx, ny );
					// build all the triangles of the fin
					build_tri_fin(returned.fanpts, tessptadj);
					
					delete returned.fanpts;
					
					// return (new edges, 0) (AFter we getnerate new edges, below)
					toret.fanpts = 0;
				}
			}
			else
			{
				// not making a fan
				tessptadj = returned.newp;
			}
			
			if (tessptadj == -1) tessptadj = add_npoint( nx, ny );
			
			nx = (NPT(ee[0], 0) + NPT(ee[1], 0)) / 2;
			ny = (NPT(ee[0], 1) + NPT(ee[1], 1)) / 2;
			tesspt = add_npoint(nx, ny);
			
			// tesselate triangles, return tesselated edge
			// due to adj_sleeve_edge_tri, ee[0] == adjee[0]
			add_newtri(exp_tri(ee[0], tesspt, tessptadj));  // tri on one side
			// quad on other side
			add_newtri(exp_tri(ee[1], tesspt, tessptadj));
			add_newtri(exp_tri(adjee[1], ee[1], tessptadj));
			
			toret.newp = tesspt;
			// return (new edges, 0)
		}
		break;
	case 3:    // joint
		{
			// we're going to need to tesselate this no matter what
			double nx, ny;
			weighted_center(	NPT(TRI(tri, 0), 0), NPT(TRI(tri, 0), 1),
				NPT(TRI(tri, 1), 0), NPT(TRI(tri, 1), 1),
				NPT(TRI(tri, 2), 0), NPT(TRI(tri, 2), 1),
				nx, ny
				);
			
			int cenpt = add_npoint( nx, ny );
			
			// for each edge e
			for (int i = 0;i < 3;i++)
			{
				int adjt = NEIGH(tri, i);
				
				exp_edge adjee;
				if ( (adjt != -1) && !shedge(tri, adjt, &adjee ) )
				{
					OT_ASSERT(0);
				}
				bool maketri = false;
				
				if ( (adjt == -1) )
				{ // can't use shedge, need to find the explict edge another way
					int countused[3];
					for (int i = 0;i < 3;i++) countused[i] = 0;
					exp_edge other;
					for (int i = 0;i < 3;i++)
					{
						if ( NEIGH(tri, i) != -1 )
						{
							OT_ASSERT(shedge(tri, NEIGH(tri, i), &other ));
							for (int j = 0;j < 2;j++)
							{
								for (int k = 0;k < 3;k++)
								{
									if (other[j] == TRI(tri, k))
									{
										countused[k]++;
									}
								}
							}
						}
					}
					int adjeecnt = 0;
					for (int i = 0;i < 3;i++)
					{
						if (countused[i] != 2)
						{
							adjee[adjeecnt++] = TRI(tri, i);
						}
					}
				}
				
				if ((adjt != -1) && !(adjee == ee))
				{ // only want to do it if it's not our parent
					// (note if adjt == -1, we act as if it's our parent, as it's faked from gen spine
					
					// 1. call tri on the other side of the edge
					returned = gen_spine_tris(adjt, adjee);
					// 2.if result has nonzero triangle fan list
					if (returned.fanpts)
					{
						// create tris from center of this tri
						// and join them to the fan list
						
						// build all the triangles of the fin
						build_tri_fin(returned.fanpts, cenpt);
						delete returned.fanpts;
					}
					else
					{
						// else if result has edges
						maketri = true;
						// the new edges now replace e
						// create new edge from center to new edges
						// create triangles to the new edges
					}
				}
				else
				{  // our parent side, need to give something to return
					// we need to tesselate the adjacent edge
					nx = (NPT(adjee[0], 0) + NPT(adjee[1], 0)) / 2;
					ny = (NPT(adjee[0], 1) + NPT(adjee[1], 1)) / 2;
					
					returned.newp = add_npoint( nx, ny, adjt == -1 );
					// return tesselated point
					toret.newp = returned.newp;
					maketri = true;
				}
				
				if (maketri)
				{
					add_newtri(exp_tri(adjee[0], returned.newp, cenpt));  // tri on one side
					add_newtri(exp_tri(adjee[1], returned.newp, cenpt));  // tri on one side
				}
			}
			// if we are the very root, ee.p1==-1
			// else we will need to return the tesselated edge we
			// were called from
			
			// for each edge e
			// 1. call tri on the other side of the edge
			// 2.if result has nonzero triangle fan list
			// create tris from center of this tri
			// and join them to the fan list
			// else if result has edges
			// the new edges now replace e
			// create new edge from center to new edges
			// create triangles to the new edges
			//
			//
				}
				break;
			}
			return toret;
}
		
void ot_inflate::generate_spine()
{
	OT_ASSERT(!spined);
	OT_ASSERT(triangulated);
	bool force_neigh = false;
	
	// find a joint first
	fjmarked = new int[out.numberoftriangles];
	memset(fjmarked, 0, out.numberoftriangles*sizeof(int));
	int root_joint = find_root_joint(0);
	if (root_joint == -1)
	{ // nobody has 3 neighbors, so pretend somebody does
		root_joint = 0;
		while ( NEIGHCNT(root_joint) == 1 )
		{
			// for each edge e
			for (int i = 0;i < 3;i++)
			{
				if ( NEIGH(root_joint, i) != -1 )
				{
					root_joint = NEIGH(root_joint, i);
					i = 3;
				}
			}
		}
		// by now root_joint should have 2 neighbors
		force_neigh = true;
	}
	
	OT_ASSERT(NEIGHCNT(root_joint) > 1);
	OT_ASSERT(root_joint != -1); // this can actually happen, so this is temporary
	
	delete[] fjmarked;
	fjmarked = NULL;
	
	// we'll need at most one new point per edge for the spine
	int max_new_points = ptct * 2 + 2 * out.numberofedges;
	npoints = new double[max_new_points];
	memset(npoints, 0, sizeof(double) * (max_new_points) );
	for (size_t i = 0;i < (ptct*2);i++)
	{
		npoints[i] = points[i];
	}
	nptct = ptct;  // start counting after our curr points
	// accumulate the sum to get the average of each spine point
	spine_sum_dist = new double[max_new_points];
	memset(spine_sum_dist, 0, max_new_points*sizeof(double));
	spine_sum_total = new int[max_new_points];
	memset(spine_sum_total, 0, max_new_points*sizeof(int));
	
	if (!force_neigh)
	{
		gen_spine_tris(root_joint, exp_edge( -1, -1) );
	}
	else
	{
		// need an extra slot for a new external point
		nptct++;
		gst_return r;
		// for each edge e
		for (int i = 0;i < 3;i++)
		{
			if ( NEIGH(root_joint, i) == -1 )
			{
				i = 3;
				r = gen_spine_tris(root_joint, exp_edge( -1, -1), 3 );
			}
		}
	}
	
	spined = true;
}
		
// half number of edges an edge is tesselated into
//#define EDGETESS 6
// now an argument (tm)

void ot_inflate::generate_3d_shape(int EDGETESS)
{
	OT_ASSERT(!made_3d);
	OT_ASSERT(triangulated);
	OT_ASSERT(spined);
	OT_ASSERT( ((EDGETESS % 2) == 0) && (EDGETESS > 0) );
	// copy each npoint to 3d points
	
	// we need nptct points,
	//plus an additional spine # of points
	//plus for each int_ext edge, add EDGETESS-2 pts (the points that occur in between the middle and the top)
	n3ptct = (nptct + (nptct - ptct) + (EDGETESS - 2) * intext_v.size() ) ;
	
	size_t max_new_doubles = 3 * n3ptct ;
	n3pts = new double[max_new_doubles];
	
	memset(n3pts, 0, sizeof(double) * (max_new_doubles) );
	
	for (size_t i = 0;i < nptct;i++)
	{
		n3pts[i*3] = npoints[i * 2];
		n3pts[i*3 + 1] = npoints[i * 2 + 1];
		if (i >= ptct)
			n3pts[i*3 + 2] = spine_sum_dist[i] / spine_sum_total[i];
	}
	// copy spine twice since there is an upper and lower
	for (size_t i = nptct;i < (nptct + (nptct - ptct));i++)
	{
		n3pts[i*3] = npoints[( (i - nptct) + ptct) * 2];
		n3pts[i*3 + 1] = npoints[( (i - nptct) + ptct) * 2 + 1];
		n3pts[i*3 + 2] = -1.0 * spine_sum_dist[i - (nptct - ptct)] / spine_sum_total[i - (nptct - ptct)];
	}
	
	// foreach int_ext edge...
	// tesselate into a constant number of points
	// each point should be moved in Z in the func sqrt(r^2 - t^2)
	for (size_t i = 0;i < intext_v.size();i++)
	{
		size_t curindex = (nptct + (nptct - ptct)) + i * (EDGETESS - 2);
		
		//second point in the edge is always the spine (internal) point (sx,sy, sz)
		double sx = N3PT(intext_v[i][1], 0);
		double sy = N3PT(intext_v[i][1], 1);
		double sz = N3PT(intext_v[i][1], 2);
		
		double ex = N3PT(intext_v[i][0], 0);
		double ey = N3PT(intext_v[i][0], 1);
		double ez = N3PT(intext_v[i][0], 2);
		
		double a = 0;
		for (int j = 0;j < (EDGETESS - 2);j++)
		{
			a += (M_PI / EDGETESS);
			if ( j == ((EDGETESS - 2) / 2) ) a += (M_PI / EDGETESS); // skip center point, since we have it already
			
			// y 'axis' is the vector from the (sx,sy,0) to (ex, ey, ez)
			//(ez==0, so we don't subtracting ez-0)
			// x 'axis' is the vector from the (sx,sy,0) to (sx, sy, sz)
			n3pts[ (j + curindex)*3 ] = sx + sin(a) * (ex - sx);
			n3pts[ (j + curindex)*3 + 1] = sy + sin(a) * (ey - sy);
			n3pts[ (j + curindex)*3 + 2] = 0 + cos(a) * sz;
		}
		
	}
	
	// foreach triangle
	for (std::list<exp_tri>::iterator i = newtri.begin();
	i != newtri.end(); ++i)
	{
		// figure out which edges are intext (go from the outside to the spine)
		int intcnt = 0;
		
		int intpts[3];
		int extpts[3];
		
		for (size_t k = 0;k < 3;k++)
		{
			if ( is_external( ((*i)[k]) ) )
			{
				extpts[k - intcnt] = (*i)[k];
			}
			else
			{
				intpts[intcnt++] = (*i)[k];
			}
		}
		OT_ASSERT( (intcnt > 0) && (intcnt < 3) );
		
		// we either have a triangle to tesselate with 2 spine points, or just 1;
		// code is specialized for either, since the total shape with both will
		// end up being either a twisted bowtie or rhombus
		if (intcnt == 2)
		{
			OT_ASSERT(intext_h.find(exp_edge(extpts[0], intpts[0])) != intext_h.end() );
			OT_ASSERT(intext_h.find(exp_edge(extpts[0], intpts[1])) != intext_h.end() );
			size_t edge0 = intext_h[exp_edge(extpts[0], intpts[0])];
			size_t edge1 = intext_h[exp_edge(extpts[0], intpts[1])];
			for (size_t j = 0;j < EDGETESS;j++)
			{
				int start = j;
				int end = j + 1;
				
				int s0, s1;
				int e0, e1;
				
				if (start == 0){ s0 = intpts[0]; s1 = intpts[1]; }
				else
				{
					if ( start != (EDGETESS / 2) )
					{
						s0 = (nptct + (nptct - ptct)) + start - 1 + edge0 * (EDGETESS - 2);
						s1 = (nptct + (nptct - ptct)) + start - 1 + edge1 * (EDGETESS - 2);
						if (start >= (EDGETESS / 2)){ s0--; s1--; }
					}
				}
				if (start == (EDGETESS - 1))
				{
					e0 = intpts[0] + (nptct - ptct);
					e1 = intpts[1] + (nptct - ptct);
				}
				else
				{
					if (end != (EDGETESS / 2))
					{
						e0 = (nptct + (nptct - ptct)) + end - 1 + edge0 * (EDGETESS - 2);
						e1 = (nptct + (nptct - ptct)) + end - 1 + edge1 * (EDGETESS - 2);
						if (end >= (EDGETESS / 2)){ e0--; e1--; }
					}
				}
				
				
				if (end == (EDGETESS / 2))
				{ // halfway point, make a single triangle
					n3tri.push_back( exp_tri(s0, s1, extpts[0]) );
				}
				else if (start == (EDGETESS / 2))
				{
					n3tri.push_back( exp_tri(extpts[0], e1, e0) );
				}
				else
				{
					n3tri.push_back( exp_tri(s0, s1, e0) );
					n3tri.push_back( exp_tri(e0, s1, e1) );
				}
			}
		}
		else if (intcnt == 1)
		{
			OT_ASSERT(intext_h.find(exp_edge(extpts[0], intpts[0])) != intext_h.end() );
			OT_ASSERT(intext_h.find(exp_edge(extpts[1], intpts[0])) != intext_h.end() );
			size_t edge0 = intext_h[exp_edge(extpts[0], intpts[0])];
			size_t edge1 = intext_h[exp_edge(extpts[1], intpts[0])];
			
			for (int j = 0;j < EDGETESS;j++)
			{
				int start = j;
				int end = j + 1;
				
				int s0, s1;
				int e0, e1;
				
				if (start == 0){ s0 = intpts[0]; s1 = intpts[1]; }
				else
				{
					if ( start != (EDGETESS / 2) )
					{
						s0 = (nptct + (nptct - ptct)) + start - 1 + edge0 * (EDGETESS - 2);
						s1 = (nptct + (nptct - ptct)) + start - 1 + edge1 * (EDGETESS - 2);
						if (start >= (EDGETESS / 2)){ s0--; s1--; }
					}
					else
					{
						s0 = extpts[0]; s1 = extpts[1];
					}
				}
				if (start == (EDGETESS - 1))
				{
					e0 = intpts[0] + (nptct - ptct);
					e1 = intpts[1] + (nptct - ptct);
				}
				else
				{
					if (end != (EDGETESS / 2))
					{
						e0 = (nptct + (nptct - ptct)) + end - 1 + edge0 * (EDGETESS - 2);
						e1 = (nptct + (nptct - ptct)) + end - 1 + edge1 * (EDGETESS - 2);
						if (end >= (EDGETESS / 2)){ e0--; e1--; }
					}
					else
					{
						e0 = extpts[0]; e1 = extpts[1];
					}
				}
				
				
				if (start == 0)
				{ // top, make a single triangle
					n3tri.push_back( exp_tri(intpts[0], e0, e1) );
				}
				else if (start == (EDGETESS - 1))
				{  // bottom, also single triangle
					n3tri.push_back( exp_tri(intpts[0] + (nptct - ptct), s1, s0 ) );
				}
				else
				{
					n3tri.push_back( exp_tri(e0, s1, s0) );
					n3tri.push_back( exp_tri(e1, s1, e0) );
				}
			}
		}
	}
	
	// done!
	
	made_3d = true;
}